///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file hamiltonian.h
///
/// \brief   header file for the pcps::Hamiltonian class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PCPS_HAMILTONIAN_HEADER
#define PCPS_HAMILTONIAN_HEADER

#include <src/pcps.h>
#include <src/operator.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   class to create, store, and access the Hamiltonian
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////

  template<class S> class Hamiltonian {

    public:

      Hamiltonian(const pcps::Input &userinp);

      typedef typename pcps::OpBase<S>::const_ptr OpBaseConstPtr; ///< pointer to base operator class

      const std::vector<OpBaseConstPtr> & ops()          const { return _ops;          } ///< access operator pointers
      const std::vector<OpBaseConstPtr> & ops_diag()     const { return _ops_diag;     } ///< access diagonal operator pointers
      const std::vector<OpBaseConstPtr> & ops_off_diag() const { return _ops_off_diag; } ///< access off-diagonal operator pointers

      const std::vector<pcps::OpRL<S>       > & opsRL()       const { return _opsRL;       } ///< access raising/lowering ops
      const std::vector<pcps::OpRLRL<S>     > & opsRLRL()     const { return _opsRLRL;     } ///< access double raising/lowering ops
      const std::vector<pcps::OpNN<S>       > & opsNN()       const { return _opsNN;       } ///< access double number ops
      const std::vector<pcps::OpN<S>        > & opsN()        const { return _opsN;        } ///< access number ops
      const std::vector<pcps::OpIdentity<S> > & opsIdentity() const { return _opsIdentity; } ///< access identity ops

      const std::vector<int> & nz_pairs() const { return _nz_pairs; } ///< access pairs with non-zero two electron integrals
      double enuc() const { return _enuc; } ///< return the nuclear repulsion energy
      bool is_qc() const { return _qc; } ///< whether or not the hamiltonian is a quantum chemistry hamiltonian

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      ///< \brief   access sites touched by imaginary time operators
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      const std::vector<std::vector<int>    > & imag_time_sites() const { return _imag_time_sites; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      ///< \brief   access imaginary time operators
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      const std::vector<std::vector<double> > & imag_time_ops() const { return _imag_time_ops; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   get a one electron integral
      /// 
      /// \param[in]  p   1st index
      /// \param[in]  q   2nd index
      ///
      /// \return  the integral \f$  h_{p q}  \f$
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S get_oei(const int p, const int q) const { return _oei[ p * _norb + q ]; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   get a modified one electron integral
      /// 
      ///    The modified one electron integrals f are definded as
      ///      \f$  f_{p q} = h_{p q} - \frac{1}{2} \sum_r  [pr|rq]  \f$
      ///    where h is the original one electron integrals and [**|**] are the two electron integrals.
      ///
      /// \param[in]  p   1st index
      /// \param[in]  q   2nd index
      ///
      /// \return  the modified integral \f$  f_{p q}  \f$
      ///
      /// \sa initialize_mei
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S get_mei(const int p, const int q) const { return _mei[ p * _norb + q ]; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   compute the compound index for accessing a two electron integral
      /// 
      /// \param[in]  p   1st index for electron coordinate 1
      /// \param[in]  q   2nd index for electron coordinate 1
      /// \param[in]  r   1st index for electron coordinate 2
      /// \param[in]  s   2nd index for electron coordinate 2
      ///
      /// \return  the compound index pqrs
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      static int tei_compound_index(const int p, const int q, const int r, const int s) {

        // declare some compound indices
        int pq, rs, pqrs;

        // compute compound index for first pair
        if (p >= q) pq = p * ( p + 1 ) / 2 + q;
        else        pq = q * ( q + 1 ) / 2 + p;

        // compute compound index for second pair
        if (r >= s) rs = r * ( r + 1 ) / 2 + s;
        else        rs = s * ( s + 1 ) / 2 + r;

        // compute overall compound index
        if (pq >= rs) pqrs = pq * ( pq + 1 ) / 2 + rs;
        else          pqrs = rs * ( rs + 1 ) / 2 + pq;

        // return the total compound index
        return pqrs;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   get a two electron integral (read only)
      /// 
      /// \param[in]  p   1st index for electron coordinate 1
      /// \param[in]  q   2nd index for electron coordinate 1
      /// \param[in]  r   1st index for electron coordinate 2
      /// \param[in]  s   2nd index for electron coordinate 2
      ///
      /// \return  the integral [pq|rs]
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S get_tei(const int p, const int q, const int r, const int s) const {
        return _tei[this->tei_compound_index(p,q,r,s)];
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief   get a two electron integral (read/write)
      /// 
      /// \param[in]  p   1st index for electron coordinate 1
      /// \param[in]  q   2nd index for electron coordinate 1
      /// \param[in]  r   1st index for electron coordinate 2
      /// \param[in]  s   2nd index for electron coordinate 2
      ///
      /// \return  the integral [pq|rs]
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S & get_tei(const int p, const int q, const int r, const int s) {
        return _tei[this->tei_compound_index(p,q,r,s)];
      }

      void write_dmrg_ints(const pcps::Input &userinp) const;

    private:

      void initialize_hubbard(const pcps::Input &userinp);
      void initialize_heisenberg(const pcps::Input &userinp);
      void initialize_from_file(const pcps::Input &userinp);
      void initialize_from_psi3(const pcps::Input &userinp);
      void initialize_from_molpro(const pcps::Input &userinp);
      void initialize_qc_from_file(const pcps::Input &userinp, const bool order_1212);
      void add_operator_from_file(const std::string &line);
      void initialize_mei(const pcps::Input &userinp);
      void initialize_nz_pairs(const pcps::Input &userinp);

      const bool _qc;                                       ///< whether or not this is a quantum chemistry hamiltonian
      const int _norb;                                      ///< number of quantum chemistry spatial orbitals
      double _enuc;                                         ///< nuclear repulsion energy
      std::vector<OpBaseConstPtr> _ops;                     ///< pointers to the hamiltonian's operators
      std::vector<OpBaseConstPtr> _ops_diag;                ///< pointers to the hamiltonian's diagonal operators
      std::vector<OpBaseConstPtr> _ops_off_diag;            ///< pointers to the hamiltonian's off-diagonal operators
      std::vector<std::vector<int>    > _imag_time_sites;   ///< sites touched by the imaginary time operators
      std::vector<std::vector<double> > _imag_time_ops;     ///< list of the imaginary time operators
      std::vector<pcps::OpRL<S>       > _opsRL;        ///< raising/lowering operators
      std::vector<pcps::OpRLRL<S>     > _opsRLRL;      ///< double raising/lowering operators
      std::vector<pcps::OpNN<S>       > _opsNN;        ///< double number operators
      std::vector<pcps::OpN<S>        > _opsN;         ///< number operators
      std::vector<pcps::OpIdentity<S> > _opsIdentity;  ///< identity operators
      std::vector<S> _oei;                             ///< one electron integrals, size _norb * _norb
      std::vector<S> _mei;                             ///< modified one electron integrals, size _norb * _norb
      std::vector<S> _tei;                             ///< two electron integrals, stored using 8-fold symmetry
      std::vector<int> _nz_pairs;                           ///< pairs of orbital indices with nonzero two electron integrals

      Hamiltonian(); ///< default constructor is disabled
      Hamiltonian(const Hamiltonian &); ///< copy constructor is disabled
      Hamiltonian &operator=(const Hamiltonian &); ///< assignment operator is disabled

  };

}

#endif
